home *** CD-ROM | disk | FTP | other *** search
Text File | 1996-12-16 | 17.3 KB | 703 lines | [TEXT/CWIE] |
- unit EDM;
-
- {Euclidian distance maps, ultimate eroded points, and watershed segmentation}
-
- interface
-
- uses
- Types, Memory, QuickDraw, QuickDrawText, Packages, Menus, Events, Fonts,
- Scrap, ToolUtils, Resources, Errors, Palettes, StandardFile, Windows,
- Controls, TextEdit, Files, Dialogs, TextUtils, Finder, MixedMode, Processes,
- globals, Utilities, Graphics, Analysis;
-
-
- procedure MakeEDM(item: integer);
-
-
- implementation
-
- type
- Image16Data = packed array[0..maxImageSize] of integer;
- Image16Ptr = ^Image16Data;
- StartsArray = array[0..255] of LongInt;
-
-
- function FindUtimatePoints(MaxEDM,item: LongInt; image2: ImageP;
- LevelStart: StartsArray; xCoordinate, yCoordinate: Image16Ptr): boolean;
- {
- Finds peaks in the EDM that contain pixels equal to or greater than
- all of their neighbors.
- }
- var
- x, y, level, rowsize, offset, i, count: LongInt;
- CoordOffset, xmax, ymax: LongInt;
- NextTick: LongInt;
- image: ImageP;
- SetPixel: boolean;
- begin
- with info^ do begin
- image := ImageP(PicbaseAddr);
- BlockMove(image, image2, PixMapSize);
- rowsize := BytesPerRow;
- xmax := PixelsPerLine - 1;
- ymax := nLines - 1;
- NextTick := TickCount + 15;
- for level := maxEDM - 1 downto 1 do begin
- repeat
- count := 0;
- for i := 0 to Histogram[level] - 1 do begin
- CoordOffset := LevelStart[level] + i;
- x:= xCoordinate^[CoordOffset];
- y:= yCoordinate^[CoordOffset];
- offset := x + y * rowsize;
- if image^[offset] <> 255 then begin
- SetPixel := false;
- if (x > 0) and (y > 0) then
- if image^[offset - rowSize - 1] > level then
- SetPixel := true;
- if y > 0 then
- if image^[offset - rowSize] > level then
- SetPixel := true;
- if (x < xmax) and (y > 0) then
- if image^[offset - rowSize + 1] > level then
- SetPixel := true;
- if x < xmax then
- if image^[offset + 1] > level then
- SetPixel := true;
- if (x < xmax) and (y < ymax) then
- if image^[offset + rowSize + 1] > level then
- SetPixel := true;
- if y < ymax then
- if image^[offset + rowSize] > level then
- SetPixel := true;
- if (x > 0) and (y < ymax) then
- if image^[offset + rowSize - 1] > level then
- SetPixel := true;
- if x > 0 then
- if image^[offset - 1] > level then
- SetPixel := true;
- if SetPixel then begin
- image^[offset] := 255;
- count := count + 1;
- end;
- end; {if}
- end; {for i}
- until count = 0;
- if TickCount > NextTick then begin
- ShowAnimatedWatch;
- if CommandPeriod then begin
- beep;
- undo;
- WhatToUndo := NothingToUndo;
- FindUtimatePoints := false;
- exit(FindUtimatePoints);
- end;
- NextTick := TickCount + 15;
- end;
- end; {for level}
- if item = WatershedItem then begin
- for i := 0 to info^.PixMapSize - 1 do
- if (image^[i] > 0) and (image^[i] < 255) then
- image2^[i] := 255;
- BlockMove(image2, image, PixMapSize);
- end else begin
- for i := 0 to info^.PixMapSize - 1 do
- if image^[i] = 255 then
- image^[i] := 0;
- end;
- FindUtimatePoints := true;
- end; {with}
- end;
-
-
- procedure DoWatershedSegmentation(MaxEDM: LongInt; image2: ImageP;
- LevelStart: StartsArray; xCoordinate, yCoordinate: Image16Ptr;
- var iterations: LongInt);
- {
- Assumes the current image contains an EDM and that the peaks in the EDM
- (the ultimate eroded points) have been set to 255. The EDM is dilated from
- these peaks, starting at the highest peak and working down. At each level,
- the dilation is constrained to pixels whose values are at that level and
- also constrained to prevent features from merging.
- }
- var
- rowSize, NextTicks: LongInt;
- level, x, y, offset, xmax, ymax: LongInt;
- count, FirstCount: LongInt;
- table: FateTable;
- image: ImageP;
-
- procedure CheckAbort;
- begin
- UpdatePicWindow;
- if CommandPeriod then begin
- beep;
- undo;
- WhatToUndo := NothingToUndo;
- exit(DoWatershedSegmentation);
- end;
- NextTicks := TickCount + 30;
- if OptionKeyDown then
- wait(30);
- end;
-
- procedure MakeFateTable;
- {Add pixel on 1st pass if bit 0 is set, on 2nd pass if bit 1 is}
- {set, on 3rd pass if bit 2 is set, and on 4th pass if bit 3 is set.}
- {E.g. '4' = add on 3rd pass, '3' = add on either 1st or 2nd pass,}
- {'f' = add on any pass. There is a routine in 'user.p' that draws}
- {all 256 neighborhoods.}
- const
- s999 = '01234567890123456789012345678901';
-
- s000 = '0044004480cc80cc0000000080cc80cc';
- s032 = '1000000080ff80ff1000000090ff90ff';
- s064 = '00000000000000000000000000000000';
- s096 = '1000000090ff90ff1000000090ff90ff';
- s128 = '2266006600ff00ff0000000000ff00ff';
- s160 = '13ff00ffffffffff33ff00ffffffffff';
- s192 = '2266006600ff00ff0000000000ff00ff';
- s224 = '33ff00ffffffffff33ff00fffffffff';
- var
- s: str255;
- c: char;
- i: LongInt;
- begin
- s := concat(s000, s032, s064, s096, s128, s160, s192, s224);
- for i := 0 to 254 do begin
- c := s[i + 1];
- if c <= '9' then
- table[i] := ord(c) - ord('0')
- else
- table[i] := 10 + ord(c) - ord('a')
- end;
- table[255] := 15; {f}
- end;
-
- procedure ProcessLevel(level, pass: LongInt);
- var
- index, x, y, i, CoordOffset, offset: LongInt;
- begin
- BlockMove(image, image2, info^.PixMapSize);
- for i := 0 to Histogram[level] - 1 do begin
- CoordOffset := LevelStart[level] + i;
- x:= xCoordinate^[CoordOffset];
- y:= yCoordinate^[CoordOffset];
- offset := x + y * rowsize;
- if image2^[offset] <> 255 then begin
- index := 0;
- if (x > 0) and (y > 0) then
- if image2^[offset - rowSize - 1] = 255 then
- index := bor(index, 1);
- if y > 0 then
- if image2^[offset - rowSize] = 255 then
- index := bor(index, 2);
- if (x < xmax) and (y > 0) then
- if image2^[offset - rowSize + 1] = 255 then
- index := bor(index, 4);
- if x < xmax then
- if image2^[offset + 1] = 255 then
- index := bor(index, 8);
- if (x < xmax) and (y < ymax) then
- if image2^[offset + rowSize + 1] = 255 then
- index := bor(index, 16);
- if y < ymax then
- if image2^[offset + rowSize] = 255 then
- index := bor(index, 32);
- if (x > 0) and (y < ymax) then
- if image2^[offset + rowSize - 1] = 255 then
- index := bor(index, 64);
- if x > 0 then
- if image2^[offset - 1] = 255 then
- index := bor(index, 128);
- case pass of
- 1: if band(table[index], 1) = 1 then begin
- image^[offset] := 255;
- count := count + 1;
- end;
- 2: if band(table[index], 2) = 2 then begin
- image^[offset] := 255;
- count := count + 1;
- end;
- 3: if band(table[index], 4) = 4 then begin
- image^[offset] := 255;
- count := count + 1;
- end;
- 4: if band(table[index], 8) = 8 then begin
- image^[offset] := 255;
- count := count + 1;
- end;
- end; {case}
- end; {if}
- end; {for i}
- end; {ProcessLevel}
-
- procedure PostProcess;
- var
- i: LongInt;
- begin
- for i := 0 to info^.PixMapSize - 1 do
- if image^[i] < 255 then
- image^[i] := 0
- end;
-
- begin
- with info^ do begin
- rowSize := BytesPerRow;
- MakeFateTable;
- image := ImageP(PicbaseAddr);
- xmax := PixelsPerLine - 1;
- ymax := nLines - 1;
- NextTicks := TickCount;
- iterations := 0;
- for level := maxEDM - 1 downto 1 do begin
- repeat
- count := 0;
- ProcessLevel(level, 1);
- ProcessLevel(level, 3);
- ProcessLevel(level, 2);
- ProcessLevel(level, 4);
- iterations := iterations + 1;
- until count = 0;
- if TickCount > NextTicks then
- CheckAbort;
- end; {for level}
- PostProcess;
- UpdatePicWindow;
- end;
- end;
-
-
- procedure SmoothEDM(image2: ImageP);
- var
- x, y, rowsize, offset, sum: LongInt;
- xmax, ymax: LongInt;
- image: ImageP;
- begin
- with info^ do begin
- image := ImageP(PicbaseAddr);
- BlockMove(image, image2, PixMapSize);
- rowsize := BytesPerRow;
- xmax := PixelsPerLine - 1;
- ymax := nLines - 1;
- for y := 0 to nLines - 1 do begin
- for x := 0 to PixelsPerLine - 1 do begin
- offset := x + y * rowsize;
- if image2^[offset] <> 1 then begin
- sum := image2^[offset] * 2;
- if (x > 0) and (y > 0) then
- sum := sum + image2^[offset - rowSize - 1];
- if y > 0 then
- sum := sum + image2^[offset - rowSize];
- if (x < xmax) and (y > 0) then
- sum := sum + image2^[offset - rowSize + 1];
- if x < xmax then
- sum := sum + image2^[offset + 1];
- if (x < xmax) and (y < ymax) then
- sum := sum + image2^[offset + rowSize + 1];
- if y < ymax then
- sum := sum + image2^[offset + rowSize];
- if (x > 0) and (y < ymax) then
- sum := sum + image2^[offset + rowSize - 1];
- if x > 0 then
- sum := sum + image2^[offset - 1];
- image^[offset] := sum div 10;
- end;
- end; {for x}
- end; {for y}
- end; {with}
- end;
-
-
- procedure MakeEDM(item: integer);
- {
- Converts a binary image into a grayscale Euclidian distance
- map (EDM). Each foreground (black) pixel in the binary image
- is assigned a value equal to its distance from the nearest
- background (white) pixel. Uses the two pass EDM algorithm
- from the "Image Pricessing Handbook" by John Russ.
- }
- const
- one = 41;
- sqrt2 = 58; {~41*sqrt(2)}
- sqrt5 = 92; {~41*sqrt(5)}
-
- var
- x, y, xmax, ymax, i, MaxEDM: LongInt;
- StartTicks, NextTicks, offset, rowsize: LongInt;
- iterations: LongInt;
- PointsOK: boolean;
- image16: Image16Ptr;
- image2: ImageP;
- str: str255;
- LevelStart: StartsArray;
- xCoordinate, yCoordinate: Image16Ptr;
-
- function ConvertToIntegers:boolean;
- var
- x, y, offset: LongInt;
- begin
- with info^ do begin
- image16 := Image16Ptr(NewPtr(rowsize * nLines * SizeOf(integer)));
- if image16 = nil then begin
- ConvertToIntegers := false;
- exit(ConvertToIntegers);
- end;
- for y := 0 to nLines - 1 do
- for x := 0 to PixelsPerLine - 1 do begin
- offset := x + y * rowsize;
- image16^[offset] := ImageP(PicBaseAddr)^[offset] * one;
- end;
- ConvertToIntegers := true;
- end;
- end;
-
- procedure ConvertToBytes;
- var
- x, y, v, round, offset: LongInt;
- begin
- round := one div 2;
- MaxEDM := 0;
- with info^ do begin
- for y := 0 to nLines - 1 do
- for x := 0 to PixelsPerLine - 1 do begin
- offset := x + y * rowsize;
- v := (image16^[offset] + round) div one;
- if v > 255 then
- v := 255;
- if v > maxEDM then
- maxEDM := v;
- ImageP(PicBaseAddr)^[offset] := v;
- end;
- end;
- end;
-
- procedure SetEdgeValue;
- var
- min, inc, v: LongInt;
- r1, r2, r3, r4, r5, offimage: LongInt;
- begin
- r1 := offset - rowsize - rowsize - 2;
- r2 := r1 + rowsize;
- r3 := r2 + rowsize;
- r4 := r3 + rowsize;
- r5 := r4 + rowsize;
- min := 32767;
-
- offimage := image16^[r3 + 2];
-
- if y < 2 then
- v := offImage + one
- else
- v := image16^[r2 + 2] + one;
- if v < min then
- min := v;
-
- if x < 2 then
- v := offImage + one
- else
- v := image16^[r3 + 1] + one;
- if v < min then
- min := v;
-
- if x > xmax then
- v := offImage + one
- else
- v := image16^[r3 + 3] + one;
- if v < min then
- min := v;
-
- if y > ymax then
- v := offImage + one
- else
- v := image16^[r4 + 2] + one;
- if v < min then
- min := v;
-
- if (x < 2) or (y < 2) then
- v := offImage + sqrt2
- else
- v := image16^[r2 + 1] + sqrt2;
- if v < min then
- min := v;
-
- if (x > xmax) or (y < 2) then
- v := offImage + sqrt2
- else
- v := image16^[r2 + 3] + sqrt2;
- if v < min then
- min := v;
-
- if (x < 2) or (y > ymax) then
- v := offImage + sqrt2
- else
- v := image16^[r4 + 1] + sqrt2;
- if v < min then
- min := v;
-
- if (x > xmax) or (y > ymax) then
- v := offImage + sqrt2
- else
- v := image16^[r4 + 3] + sqrt2;
- if v < min then
- min := v;
-
- if (x < 2) or (y < 2) then
- v := offImage + sqrt5
- else
- v := image16^[r1 + 1] + sqrt5;
- if v < min then
- min := v;
-
- if (x > xmax) or (y < 2) then
- v := offImage + sqrt5
- else
- v := image16^[r1 + 3] + sqrt5;
- if v < min then
- min := v;
-
- if (x > xmax) or (y < 2) then
- v := offImage + sqrt5
- else
- v := image16^[r2 + 4] + sqrt5;
- if v < min then
- min := v;
-
- if (x > xmax) or (y > ymax) then
- v := offImage + sqrt5
- else
- v := image16^[r4 + 4] + sqrt5;
- if v < min then
- min := v;
-
- if (x > xmax) or (y > ymax) then
- v := offImage + sqrt5
- else
- v := image16^[r5 + 3] + sqrt5;
- if v < min then
- min := v;
-
- if (x < 2) or (y > ymax) then
- v := offImage + sqrt5
- else
- v := image16^[r5 + 1] + sqrt5;
- if v < min then
- min := v;
-
- if (x < 2) or (y > ymax) then
- v := offImage + sqrt5
- else
- v := image16^[r4] + sqrt5;
- if v < min then
- min := v;
-
- if (x < 2) or (y < 2) then
- v := offImage + sqrt5
- else
- v := image16^[r2] + sqrt5;
- if v < min then
- min := v;
-
- image16^[offset] := min;
- end; {SetEdgeValue}
-
- procedure SetValue;
- var
- min, inc, v: LongInt;
- r1, r2, r3, r4, r5: LongInt;
- begin
- r1 := offset - rowsize - rowsize - 2;
- r2 := r1 + rowsize;
- r3 := r2 + rowsize;
- r4 := r3 + rowsize;
- r5 := r4 + rowsize;
- min := 32767;
-
- v := image16^[r2 + 2] + one;
- if v < min then
- min := v;
- v := image16^[r3 + 1] + one;
- if v < min then
- min := v;
- v := image16^[r3 + 3] + one;
- if v < min then
- min := v;
- v := image16^[r4 + 2] + one;
- if v < min then
- min := v;
-
- v := image16^[r2 + 1] + sqrt2;
- if v < min then
- min := v;
- v := image16^[r2 + 3] + sqrt2;
- if v < min then
- min := v;
- v := image16^[r4 + 1] + sqrt2;
- if v < min then
- min := v;
- v := image16^[r4 + 3] + sqrt2;
- if v < min then
- min := v;
-
- v := image16^[r1 + 1] + sqrt5;
- if v < min then
- min := v;
- v := image16^[r1 + 3] + sqrt5;
- if v < min then
- min := v;
- v := image16^[r2 + 4] + sqrt5;
- if v < min then
- min := v;
- v := image16^[r4 + 4] + sqrt5;
- if v < min then
- min := v;
- v := image16^[r5 + 3] + sqrt5;
- if v < min then
- min := v;
- v := image16^[r5 + 1] + sqrt5;
- if v < min then
- min := v;
- v := image16^[r4] + sqrt5;
- if v < min then
- min := v;
- v := image16^[r2] + sqrt5;
- if v < min then
- min := v;
-
- image16^[offset] := min;
- end; {SetValue}
-
-
- procedure CheckAbort;
- begin
- ShowAnimatedWatch;
- if CommandPeriod then begin
- beep;
- DisposePtr(ptr(image16));
- exit(MakeEDM);
- end;
- NextTicks := TickCount + 15;
- end;
-
-
- procedure MakeCoordinateArrays;
- {Generates the XY coordinate arrays that allow pixels at each level to
- be accesed directly without searching through the entire image. This
- method, suggested by Stein Roervik (steinr@kjemi.unit.no), greatly
- speeds up the watershed segmentation routine.}
- var
- i, x, y, rowsize, offset, v, ArraySize: LongInt;
- LevelOffset: array[0..255] of LongInt;
- image: ImageP;
- begin
- with info^ do begin
- RoiRect := PicRect;
- GetRectHistogram;
- ArraySize := 0;
- for i := 1 to maxEDM - 1 do
- ArraySize := ArraySize + histogram[i];
- xCoordinate := Image16Ptr(NewPtr(ArraySize * SizeOf(integer)));
- yCoordinate := Image16Ptr(NewPtr(ArraySize * SizeOf(integer)));
- if (xCoordinate = nil) or (yCoordinate = nil) then begin
- if xCoordinate <> nil then
- DisposePtr(ptr(xCoordinate));
- DisposePtr(ptr(image2));
- PutError('Out of memory');
- undo;
- WhatToUndo := NothingToUndo;
- exit(MakeEDM);
- end;
- image := ImageP(PicbaseAddr);
- offset := 0;
- for i := 0 to 255 do begin
- LevelStart[i] := offset;
- if (i >= 1) and (i <= (MaxEDM - 1)) then
- offset := offset + histogram[i];
- end;
- for I := 0 to 255 do
- LevelOffset[i] := 0;
- rowsize := BytesPerRow;
- for y := 0 to nLines - 1 do begin
- for x := 0 to PixelsPerLine - 1 do begin
- v := image^[x + y * rowsize];
- if (v >= 1) and (v <= (MaxEDM - 1)) then begin
- offset := LevelStart[v] + LevelOffset[v];
- xCoordinate^[offset] := x;
- yCoordinate^[offset] := y;
- LevelOffset[v] := LevelOffset[v] + 1;
- end;
- end; {for x}
- end; {for y}
- end; {with}
- end;
-
- begin
- with info^ do begin
- ShowWatch;
- KillRoi;
- rowsize := BytesPerRow;
- xmax := PixelsPerLine - 3;
- ymax := nLines - 3;
- StartTicks := TickCount;
- NextTicks := TickCount;
- SetupUndo;
- WhatToUndo := UndoFilter;
- if not ConvertToIntegers then begin
- AbortMacro;
- PutError('Out of Memory');
- exit(MakeEDM);
- end;
- for y := 0 to nLines - 1 do begin
- for x := 0 to PixelsPerLine - 1 do begin
- offset := x + y * rowsize;
- if image16^[offset] > 0.0 then begin
- if (x < 2) or (x > xmax) or (y < 2) or (y > ymax) then
- SetEdgeValue
- else
- SetValue;
- end;
- end;
- if TickCount > NextTicks then
- CheckAbort;
- end;
- for y := nLines - 1 downto 0 do begin
- for x := PixelsPerLine - 1 downto 0 do begin
- offset := x + y * rowsize;
- if image16^[offset] > 0.0 then begin
- if (x < 2) or (x > xmax) or (y < 2) or (y > ymax) then
- SetEdgeValue
- else
- SetValue;
- end;
- end;
- if TickCount > NextTicks then
- CheckAbort;
- end;
- ConvertToBytes;
- DisposePtr(ptr(image16));
- if (item = UltimateItem) or (item = WatershedItem) then begin
- image2 := ImageP(NewPtr(PixMapSize));
- if image2 = nil then
- exit(MakeEDM);
- if ((item = WatershedItem) and not OptionKeyWasDown) or ((item = UltimateItem) and OptionKeyWasDown) then
- SmoothEDM(image2);
- MakeCoordinateArrays;
- PointsOK := FindUtimatePoints(MaxEDM, item, image2, LevelStart, xCoordinate, yCoordinate);
- end;
- if (item = WatershedItem) and PointsOK then begin
- DoWatershedSegmentation(MaxEDM, image2, LevelStart, xCoordinate, yCoordinate, iterations);
- str := StringOf(MaxEDM - 1:1, ' levels', crStr, iterations/(MaxEDM - 1):1:2, ' iterations/level');
- end else
- str := '';
- if (item = UltimateItem) or (item = WatershedItem) then begin
- DisposePtr(ptr(xCoordinate));
- DisposePtr(ptr(yCoordinate));
- DisposePtr(ptr(image2));
- end;
- ShowTime(StartTicks, PicRect, str);
- changes := true;
- UpdatePicWindow;
- end; {with}
- end;
-
- end.